        !! module utilities_module
        !! for doing little common things
      module utilities_module
        contains
            subroutine mypause()
            print *, " Press <enter> to continue"
            read(*,*)
            end subroutine mypause
      end module utilities_module

! to compile:
! gfortran -o umfpackInput.exe -ffixed-line-length-132 umfpackInput.f umfpack_2.0_prompt.f
! or transform to double precision
! gfortran -o umfpackInput.exe -freal-4-real-8 -ffixed-line-length-132 umfpackInput.f umfpack_2.0_prompt.f
!
! to run:
! ./umfpackInput.exe

      program umfpackInput
         use utilities_module
c*********************************************************************72
c
c umfTest is the main program for UMFPACK_2.0_PRB.
c
c    umfTest tests the UMFPACK.f
c
c    umfpack demo program.
c
c
c  Modified:
c
c    11 April 2022
c
      Parameter(NCHANL=7) ! for output
      Parameter(NZ=2500000)
      Parameter(NYDIM=900000)
      Parameter(LICN=5*NZ,LIRN=2*NZ)
      integer nmax, nemax, lvalue, lindex
      parameter (nmax=NYDIM, nemax=NZ, lvalue=LICN, lindex=LIRN)
      integer keep (20), index (lindex), info (40),
     &  i, icntl (20), n, ne, ai (2*nemax)
      real b (nmax), x (nmax), w (4*nmax), value (lvalue),
     & cntl (10), rinfo (20), ax (nemax)
      real Ajac(lvalue), Asave(lvalue) ! for small tests, otherwise use large lvalue
      integer  IRS(lindex), ICS(lindex) ! for small tests, otherwise use large lindex
      character(len=80) :: FileIn, string
      logical F01

      open(NCHANL,file='sparsmat.tst',status='unknown')
       FileIn='testMatrix.out'
!       FileIn='testIMSLmatrixTHR0.out'
      FileIn='YudintstMatrix.out'
      open(327,file=FileIn(1:len_trim(FileIn)),status='old')

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'umfTest:'
      write ( *, '(a)' ) 'Test umfpackTest.f.'
      write (NCHANL, '(a)' ) 'Test umfpackTest.f.'
c
c  read input matrix and define right-hand side.
c
c  keep a copy of the triplet form in ai and ax.
c

      N=0
      nza=0


!: read the input matrix *

!-!    read(327,'(1p,2e15.3,i10,l10)') PIVOT, THR, MAXIT, F01
      read(327,*) PIVOT, THR, MAXIT, F01
      read(327,'(a)')string ! '#    A                         IRN                       ICN'

      idd=0
      N=0
      do
        idd=idd+1
        read(327,*,end=11) Asave(idd), IRS(idd), ICS(idd)
        Ajac(idd)=Asave(idd)
        N=max(N,IRS(idd))
      enddo
11    NZA=idd-1
      write(*,'(a,1p,g25.14,2i10)')'last line:', Asave(NZA), IRS(NZA), ICS(NZA)
      write(*,*)'all NZA =', NZA,'  lines are read'
      write(*,*)'N=',N
      write(NCHANL,*)'N=',N
!       call mypause

      ne=NZA

      do i = 1, ne
        index(i) = IRS(i)
        index(ne+i) = ICS(i)
        value(i) = Asave(i)
        ai(i) = index (i)
        ai(ne+i) = index (ne+i)
        ax(i) = value (i)
      enddo

      write(*,*)'matrix size n=',n
      write(*,*)'number of nonzero elements ne=',ne
      write(*,*) '#    A                   IRN         ICN',
     &' in triplet COO format'
      write(NCHANL,*)'matrix size n=',n
      write(NCHANL,*)'number of nonzero elements ne=',ne
      write(NCHANL,*) '#    A                   IRN         ICN',
     &' in triplet COO format'
      do i = 1, ne
        write(*,*) value (i), index(i), index(ne+i)
        write(NCHANL,*) value (i), index(i), index(ne+i)
      end do
!       call mypause
      write(*,*) '#    RHS b'
      do i = 1, n
        b(i)=0.
      end do
      do K=1,NZA ! ne
        b(IRS(K))=b(IRS(K))+AJAC(K) ! !- * 1.
      enddo
      do i = 1, n
        write(*,*) b(i)
      end do

!       call mypause

c
c initialize controls, and change default printing control.  note that
c this change from the default should only be used for test cases.  it
c can generate a lot of output for large matrices.
c
      call ums2in (icntl, cntl, keep)
      icntl (3) = 4
c
c  Factorize a, and print the factors.  input matrix is not preserved.
c
      job=0 ! input matrix is not preserved.
      call ums2fa (n, ne, job, .false., lvalue, lindex, value, index,
     &               keep, cntl, icntl, info, rinfo)
      if (info (1) .lt. 0) then
        stop
      end if
c
c  Reset default printing control (ums2in could be called instead)
c
      icntl (3) = 2
c
c  Solve ax = b without iterative refinement, and print solution.
c
      icntl (8) = 0
      call ums2so (n, 0, .false., lvalue, lindex, value, index,
     &               keep, b, x, w, cntl, icntl, info, rinfo)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  begin Solution of A * x = b:'
      write ( *, '(a)' ) ' '
      write (*, 30) (x (i), i = 1, n)
      write ( *, '(a)' ) '  end Solution of A * x = b'
      write (NCHANL, '(a)' ) ' '
      write (NCHANL, '(a)' ) '  begin Solution of A * x = b:'
      write (NCHANL, '(a)' ) ' '
      write (NCHANL, 30) (x (i), i = 1, n)
      write (NCHANL, '(a)' ) '  end Solution of A * x = b'
      ERREST=0.
      do I=1,n;
         ERREST=ERREST+ABS(1.D0-X(I))
       enddo;
      write(*,'(a)' )'for unity solutions only:'
      write(*,'(a,1p,e10.3/a,e10.3)' )
     &' EXACT TOTAL ERROR=',ERREST,' average error=',ERREST/float(n)
      write(NCHANL,'(a)')'for unity solutions only:'
      write(NCHANL,'(a,1p,e10.3/a,e10.3)')
     &' EXACT TOTAL ERROR=',ERREST,' average error=',ERREST/float(n)


!       call mypause

      if (info (1) .lt. 0) then
        stop
      end if
c
c
c  Solve ax = b with iterative refinement, and print solution.
c  This is much more accurate.

      ne=NZA

      do i = 1, ne
        index(i) = IRS(i)
        index(ne+i) = ICS(i)
        value(i) = Asave(i)
        ai(i) = index (i)
        ai(ne+i) = index (ne+i)
        ax(i) = value (i)
      enddo

      write(*,*)'matrix size n=',n
      write(*,*)'number of nonzero elements ne=',ne
      write(*,*) '#    A                   IRN         ICN',
     &' in triplet COO format'
      write(NCHANL,*)'matrix size n=',n
      write(NCHANL,*)'number of nonzero elements ne=',ne
      write(NCHANL,*) '#    A                   IRN         ICN',
     &' in triplet COO format'
      do i = 1, ne
        write(*,*) value (i), index(i), index(ne+i)
        write(NCHANL,*) value (i), index(i), index(ne+i)
      end do
!       call mypause
      write(*,*) '#    RHS b'
      do i = 1, n
        b(i)=0.
      end do
      do K=1,NZA ! ne
        b(IRS(K))=b(IRS(K))+AJAC(K) ! !- * 1.
      enddo
      do i = 1, n
        write(*,*) b(i)
      end do

      job=1 ! input matrix is preserved.
      call ums2fa (n, ne, job, .false., lvalue, lindex, value, index,
     &               keep, cntl, icntl, info, rinfo)


      icntl (8) = 10
      call ums2so (n, 0, .false., lvalue, lindex,  value, index,
     &               keep, b, x, w, cntl, icntl, info, rinfo)
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Solution of A * x = b'
      write ( *, '(a)' ) '  via iterative refinement'
      write ( *, '(a)' ) ' '
      write (*, 30) (x (i), i = 1, n)
      write (*,*)'end solution of A via iterative refinement'
      write (NCHANL, '(a)' ) ' '
      write (NCHANL, '(a)' ) '  begin Solution of A * x = b:'
      write (NCHANL, '(a)' ) '  via iterative refinement'
      write (NCHANL, '(a)' ) ' '
      write (NCHANL, 30) (x (i), i = 1, n)
      write (NCHANL,*)'end solution of A via iterative refinement'
      ERREST=0.
      do I=1,n;
         ERREST=ERREST+ABS(1.D0-X(I))
       enddo;
      write(*,'(a)' )'for unity solutions only:'
      write(*,'(a,1p,e10.3/a,e10.3)' )
     &' EXACT TOTAL ERROR=',ERREST,' average error=',ERREST/float(n)
      write(NCHANL,'(a)')'for unity solutions only:'
      write(NCHANL,'(a,1p,e10.3/a,e10.3)')
     &' EXACT TOTAL ERROR=',ERREST,' average error=',ERREST/float(n)


      if (info (1) .lt. 0) then
        stop
      end if

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) ' umfTest:'
      write ( *, '(a)' ) '  Normal end of execution'

      stop
30    format ('begin solution: ', 10(/,f20.16))
      end program umfpackInput
